Poravnanja - 2. deo¶
Algoritmi za izračunavanje poravnanja sekvenci, predstavljeni u prošlom delu, mogu se dodatno unaprediti. Jedno od unapređenja ogleda se u boljoj vremenskoj i prostornoj složenosti algoritama. Svi, do sada, predstavljeni algoritmi imaju kvadratnu složenost $O(nm)$, gde je $n$ dužina prve a $m$ dužina druge sekvence.
Jedan pokušaj optimizacije može biti smanjenje prostorne složenosti, pamćenjem samo rezultata iz prethodne iteracije. Ukoliko su tabele dinamičkog programiranja bile popunjavane po redovima, dovoljno je pamtiti tekući i prethodni red tabele, u suprotnom, pamtiće se tekuća i prethodna kolona. Ovim pristupom, na prvi pogled, gubi se mogućnost rekonstrukcije poravnanja - računa se samo konačni skor poravnanja. Vremenska složenost algoritma i dalje ostaje kvadratna.
Pametnim pristupom izračunavanju, moguće je modifikovati algoritme tako da njihova prostorna složenost postane linearna, $O(min\{m,n\})$ (vremenska i dalje kvadratna) a da pri tom rekonstrukcija poravnanja bude moguća.
Primer efikasnijeg pristupa biće prikazan pomoću Hirshberg-ovog algoritma za optimizaciju Needleman-Wunsch algoritma
Dodatni materijali: Hirshberg's algorithm
Originalni Needleman-Wunsch Algoritam
MATCH_SCORE = 2
MISSMATCH_PENALTY = -1
GAP_PENALTY = -2
def score(c1, c2):
if c1 == c2:
return MATCH_SCORE
else:
return MISSMATCH_PENALTY
def needleman_wunsch(v, w):
n = len(v) + 1
m = len(w) + 1
s = [[0 for j in range(m)] for i in range(n)]
backtrack = [[None for j in range(m)] for i in range(n)]
for i in range(1, n):
s[i][0] = s[i - 1][0] + GAP_PENALTY
backtrack[i][0] = (i - 1, 0)
for j in range(1, m):
s[0][j] = s[0][j - 1] + GAP_PENALTY
backtrack[0][j] = (0, j - 1)
for i in range(1, n):
for j in range(1, m):
from_top = s[i - 1][j] + GAP_PENALTY
from_left = s[i][j - 1] + GAP_PENALTY
from_diagonal = s[i - 1][j - 1] + score(v[i - 1], w[j - 1])
s[i][j] = max(from_top, from_left, from_diagonal)
if s[i][j] == from_top:
backtrack[i][j] = (i - 1, j)
elif s[i][j] == from_left:
backtrack[i][j] = (i, j - 1)
else:
backtrack[i][j] = (i - 1, j - 1)
i = n - 1
j = m - 1
v_align = ''
w_align = ''
while backtrack[i][j] != None:
if backtrack[i][j] == (i - 1, j):
v_align = v[i - 1] + v_align
w_align = '-' + w_align
elif backtrack[i][j] == (i, j - 1):
v_align = '-' + v_align
w_align = w[j - 1] + w_align
else:
v_align = v[i - 1] + v_align
w_align = w[j - 1] + w_align
(i, j) = backtrack[i][j]
return v_align, w_align
Za početak, potrebno je implementirati Needleman-Wunsch algoritam tako da skor računa u linearnoj prostornoj složenosti. Vrednosti pojedinačnih poravnanja biće fiksirane dok će povratna vrednost algoritma biti poslednji red tabele dinamičkog programiranja.
import copy
def needleman_wunsch_last_line(v, w):
n = len(v) + 1
m = len(w) + 1
s = [[0 for j in range(m)] for i in range(2)]
for j in range(1, m):
s[0][j] = s[0][j - 1] + GAP_PENALTY
for i in range(1, n):
s[1][0] = i * GAP_PENALTY
for j in range(1, m):
from_top = s[0][j] + GAP_PENALTY
from_left = s[1][j - 1] + GAP_PENALTY
from_diagonal = s[0][j - 1] + score(v[i - 1], w[j - 1])
s[1][j] = max(from_top, from_left, from_diagonal)
s[0] = copy.copy(s[1])
return s[1]
Primer¶
Neka je potrebno izračunati skor globalnog poravnanja sekvenci V = AGTACGCA i W = TATGC
Tabela koju konstruiše Needleman-Wunsch algoritam je:
T A T G C
0 -2 -4 -6 -8 -10
A -2 -1 0 -2 -4 -6
G -4 -3 -2 -1 0 -2
T -6 -2 -4 0 -2 -1
A -8 -4 0 -2 -1 -3
C -10 -6 -2 -1 -3 1
G -12 -8 -4 -3 1 -1
C -14 -10 -6 -5 -1 3
A -16 -12 -8 -7 -3 1Dok je rezultat poslednjeg reda tabele, dobijen kao rezultat prethodno navedenog algoritma linearne prostorne složenosti:
v = 'AGTACGCA'
w = 'TATGC'
needleman_wunsch_last_line(v,w)
Način na koji Hirshberg-ov algoritam omogućava rekonstrukciju poravnanja, u linearnoj prostornoj složenosti, zasniva se na podeli-pa-vladaj strategiji konstrukcije algoritama.
Za početak, umesto izračunavanja cele tabele, računa se poslednji red (ili kolona, u zavisnosti od pristupa implementaciji) polovine tabele dinamičkog programiranja.
Primer
AGTACGCA = AGTA + CGCA
T A T G C
0 -2 -4 -6 -8 -10
A -2 -1 0 -2 -4 -6
G -4 -3 -2 -1 0 -2
T -6 -2 -4 0 -2 -1
A -8 -4 0 -2 -1 -3Nakon izračunavanja poslednjeg reda prve (gornje) polovine tabele, rezultat je:
[-8, -4, 0, -2, -1, -3]Druga polovina tabele se takođe može rekonstruisati u linearnoj složenosti krećući se u suprotnom smeru (od donjeg desnog ugla ka gornjem levom). Kako bi se to ostvarilo, potrebno je izračunati poravnanja obrnutih sekvenci.
Primer
AGTACGCA = AGTA + CGCAUmesto poravnanja CGCA sa sekvencom TATGC, poravnaće se obrnuta sekvenca ACGC sa obrnutom sekvencom CGTAT. Rezultat ovog poravnanja je tabela:
C G T A T
0 -2 -4 -6 -8 -10
A -2 -1 -3 -5 -4 -6
C -4 0 -2 -4 -6 -5
G -6 -2 2 0 -2 -4
C -8 -4 0 1 -1 -3Poslednji red ove tabele, obrnut kako bi se ispravilo prethodno obrtanje sekvenci za potrebe izračunavanja, je:
[-8, -4, 0, 1, -1, -3] -> [-3, -1, 1, 0, -4, -8]Sada je potrebno sagledati dobijene rezultate iz prethodnih koraka. Poslednji redovi prve i druge polovine izračunavanja su:
[-8, -4, 0, -2, -1, -3]
[-3, -1, 1, 0, -4, -8]Sumiranjem izračunatih redova dobija se:
[-8, -4, 0, -2, -1, -3]
[-3, -1, 1, 0, -4, -8] +
--------------------------------
[-11, -5, 1, -2, -5, -11]Maksimum rezultujućeg niza predstavlja poziciju koja maksimizira poravnanja obe polovine, pa će i globalno poravnanje polaznih sekvenci sigurno prolaziti kroz tu poziciju.
[-11, -5, *1, -2, -5, -11]Neka je pozicija maksimuma označena indeksom j. Poravnanje se rekurzivno nastavlja na nivou polovina $V[0:\frac{n}{2}]$ sa $W[0:j]$ i $V[\frac{n}{2}:]$ sa $W[j:]$
Dolaskom do slučaja kada je deljenjem jedna od sekvenci postala prazna, izvršava se bazni slučaj rekurzije, koji podrazumeva da se preostali karakteri neprazne sekvence poravnavaju sa prazninom. U slučaju kada jedna od niski ima dužinu $1$, primenjuje se originalni Needleman-Wunsch.
Konačni rezultat poravnanja dobija se spajanjem rezultata iz listova stabla rekurzije.
def hirschberg(v, w):
n = len(v)
m = len(w)
if n == 0:
v_align = '-' * m
w_align = w
return v_align, w_align
if m == 0:
v_align = v
w_align = '-' * n
return v_align, w_align
if n == 1 or m == 1:
(v_align, w_align) = needleman_wunsch(v, w)
return needleman_wunsch(v, w)
mid_i = n // 2
last_line_L = needleman_wunsch_last_line(v[:mid_i], w)
last_line_R = needleman_wunsch_last_line(v[mid_i::][::-1], w[::-1])
last_line_R.reverse()
max_j = None
max_value = float('-inf')
for j in range(m):
current_value = last_line_L[j] + last_line_R[j]
if current_value > max_value:
max_j = j
max_value = current_value
(v_align_L, w_align_L) = hirschberg(v[:mid_i], w[:max_j])
(v_align_R, w_align_R) = hirschberg(v[mid_i:], w[max_j:])
v_align = v_align_L + v_align_R
w_align = w_align_L + w_align_R
return v_align, w_align
v_align, w_align = hirschberg(v, w)
print(v_align)
print(w_align)
Provera u odnosu na originalni Needleman-Wunsch algoritam:
needleman_wunsch(v, w)
U praksi se navedena poravnanja koriste za poravnavanje kraćih sekvenci, dok se sa poravnanja čitavih genoma koriste heuristički pristupi poravnanja (npr. BLAST)
Ponekad kada se pojavi više uzastopnih praznina u poravnanju nema razloga da cena svake praznine nosi istu cenu. Drugim rečima, niz od 100 praznina nije potrebno kazniti 10 puta oštrije od niza od 10 praznina. Iz tog razloga, primenjuju se afina poravnanja koja jednom cenom kažnjavaju otvaranje praznine, dok je cena svake sledeće praznine u nizu dosta manja u odnosu na početnu. Time se obezbeđuje manja fragmentacija poravnanja (teži se poravnanju kod koga su praznine su grupisane).
Poravnanje sa afinim računanjem cene praznine¶
Umesto korišćenja jedne tabele dinamičkog poravnanja, za potrebe afinog poravnanja koristiće se tri tabele. Intuitivno, svaka od tabela odgovara kretanjem određenom operacijom (insercija, delecija, uparivanje karaktera). "Donja" tabela odgovara poravnavanju dodavanjem praznina u prvoj sekvenci, "srednja" tabela odgovara uparivanju karaktera i "gornja" poravnavanju dodavanjem praznina u drugoj sekvenci. Kroz "donju" tabelu dozvoljeno je kretanje (prenos rezultata) samo sa polja iznad (insercija), kroz "srednju" tabelu dozvoljeno je kretanje sa prethodne dijagonalne pozicije dok se kroz "gornju" tabelu može kretati samo sa leve strane (delecija). Skok sa srednje tabele na bilo donju ili gornju nosi cenu $\sigma$ dok svaki korak kroz gornju ili donju tabelu košta dodatno $\epsilon$. Cena kretanje kroz srednju tabelu zavisi od uparivanja odgovarajućih karaktera, ali se razmatra i moguć dolazak na srednju tabelu sa neke od preostale dve. U svakom koraku bira se rezultat koji maksimizuje rezultat.
Inicijalno su prvi red gornje tabele
Dodatni materijali: Penalizing Insertions and Deletions in Sequence Alignment
S = -3 # Cena početne praznine
E = -1 # Cena svake sledeće praznine
MATCH_SCORE = 2 # Cena uspešnog uparivanja
MISSMATCH_PENALTY = -0.5 # Cena pogrešnog uparivanja
def score(c1, c2):
if c1 == c2:
return MATCH_SCORE
else:
return MISSMATCH_PENALTY
def affine_gap_alignment(v, w):
n = len(v) + 1
m = len(w) + 1
s_upper = [[0 for j in range(m)] for i in range(n)]
s_middle = [[0 for j in range(m)] for i in range(n)]
s_lower = [[0 for j in range(m)] for i in range(n)]
backtrack = [[None for j in range(m)] for i in range(n)]
for i in range(1, n):
backtrack[i][0] = (i - 1, 0)
for j in range(1, m):
backtrack[0][j] = (0, j - 1)
for i in range(1, n):
for j in range(1, m):
s_lower[i][j] = max(
s_lower[i - 1][j] + E,
s_middle[i - 1][j] + S
)
s_upper[i][j] = max(
s_upper[i][j - 1] + E,
s_middle[i][j - 1] + S
)
s_middle[i][j] = max(
s_lower[i][j],
s_upper[i][j],
s_middle[i - 1][j - 1] + score(v[i - 1], w[j - 1])
)
if s_middle[i][j] == s_middle[i - 1][j - 1] + score(v[i - 1], w[j - 1]):
backtrack[i][j] = (i - 1, j - 1)
elif s_middle[i][j] == s_lower[i][j]:
backtrack[i][j] = (i - 1, j)
else:
backtrack[i][j] = (i, j - 1)
i = n - 1
j = m - 1
v_align = ''
w_align = ''
while backtrack[i][j] != None:
if backtrack[i][j] == (i - 1, j):
v_align = v[i - 1] + v_align
w_align = '-' + w_align
elif backtrack[i][j] == (i, j - 1):
v_align = '-' + v_align
w_align = w[j - 1] + w_align
else:
v_align = v[i - 1] + v_align
w_align = w[j - 1] + w_align
(i, j) = backtrack[i][j]
return v_align, w_align
v = 'TTCTA'
w = 'AGCCCTACA'
# v = 'AGTACGCA'
# w = 'TATGC'
v_align, w_align = affine_gap_alignment(v, w)
print(v_align)
print(w_align)
v = 'TTCTA'
w = 'AGCCCTACA'
needleman_wunsch(v, w)